#pragma once

#include "base/EulerSystem.h"

struct IsoVortex : public InitialCondition
{
	struct ExecFunctor
	{
		VTKM_EXEC_CONT ExecFunctor() :
			_value(0.12)
		{
		}

		VTKM_EXEC_CONT ExecFunctor(FloatDefault value) :
			_value(value)
		{
		}

		VTKM_EXEC void SetUp()
		{

		}

		VTKM_EXEC
			vtkm::Vec<FloatDefault, 5> valueExact(FloatDefault t, const Vec3f& p)
		{
			Vec3f _vel{ 1, 1, 0 };
			FloatDefault x = p[0] - _vel[0] * t;
			FloatDefault y = p[1] - _vel[1] * t;
			FloatDefault z = p[2] - _vel[2] * t;
			FloatDefault gam = 1.4, gamm1 = gam - 1, epi = 5.;
			FloatDefault xb, yb, r2;

			FloatDefault PI = vtkm::Pi();
			xb = x + 5;
			yb = y + 5;
			r2 = xb * xb + yb * yb;

			FloatDefault u = _vel[0] + epi / (2. * PI) * exp(0.5 * (1. - r2)) * (-yb);
			FloatDefault v = _vel[1] + epi / (2. * PI) * exp(0.5 * (1. - r2)) * xb;
			FloatDefault T = 1. - gamm1 * epi * epi / (8 * gam * PI * PI) * exp(1. - r2);
			FloatDefault density = pow(T, 1 / gamm1);
			FloatDefault pre = pow(density, gam);

			vtkm::Vec<FloatDefault, 5> value;
			value[0] = density;
			value[1] = density * u;
			value[2] = density * v;
			value[3] = density * 0;
			value[4] = pre / gamm1 + 0.5 * density * (u * u + v * v);

			return value;
		}


		VTKM_EXEC vtkm::Vec<FloatDefault, 5> operator()(const Vec3f& point)
		{
			return valueExact(0, point);
		}

		FloatDefault _value;
	};

	using ExecObjectType = ExecFunctor;

private:
	vtkm::FloatDefault _value = 100;
};

class IsoVortexProblem : public EulerSystem
{
	using Superclass = EulerSystem;
public:
	IsoVortexProblem(Mesh& mesh) :
		Superclass(mesh)
	{}

	virtual void Projection() override
	{
		InitialCondition::Execute(IsoVortex{}, _mesh.GetCellSet(), _mesh.GetCoordinate().GetDataAsMultiplexer(), _group_variable);
	}
	~IsoVortexProblem() {};

private:

};

